% Pf, future prices of individuals in each 
% costs depend on age given the state

function Ptm = fixedpointTmx_v6b(CAPH,MU,rho,Pf,Ptm,tolFP,maxreps)

%number of future periods
tf = size(Pf,3);
T = size(CAPH,3)+1;
% this includes death as state
Nstates = size(CAPH,1);
if isempty(Ptm)
	Ptm = zeros(Nstates,1);
end
inP = Ptm + 2*tolFP;
reps = 1;

% discount from now until tf
t0 = T -tf;
% future costs
muf = MU(:,:,t0+1:T);
% current costs
mu0 = MU(:,:,t0);
while (max(abs(inP-Ptm)>tolFP)==1 && reps<maxreps)    
    inP = Ptm;
    % Pf future of P, is an N x N x (T-t) matrix
    % dummy of staying conditional on getting 
    % to the period in question
    stay  = repmat(permute(Pf,[2,1,3]),Nstates,1)>=repmat(inP,1,Nstates,tf);
    % death is a no stay state
    stay(:,Nstates,:) = 0;
    % probability of getting to each state in a 
    % given period -- transit from periods where i stay 
    capH=CAPH(:,:,t0);
    pt = zeros(Nstates,Nstates,tf);
    pt(:,:,1) = capH;
    E = zeros(Nstates,Nstates,tf);
    
    E(:,:,1) = rho*pt(:,:,1).*stay(:,:,1).*repmat(permute(muf(:,:,1),[2,1,3]),Nstates,1,1);

    D = zeros(Nstates,Nstates,tf);
    D(:,:,1)= rho*pt(:,:,1).*stay(:,:,1);
    for t=2:tf
        capH=CAPH(:,:,t0+t-1);
        pt(:,:,t) = pt(:,:,t-1).*stay(:,:,t-1)*capH;
        E(:,:,t) = rho^t*pt(:,:,t).*stay(:,:,t).*repmat(permute(muf(:,:,t),[2,1,3]),Nstates,1,1);
        D(:,:,t)= rho^t*pt(:,:,t).*stay(:,:,t);
    end

    % numerator
    E = sum(sum(E,2),3);
    num = mu0+E;
    % denominator
    den = 1 + sum(sum(D,2),3);
    Ptm = num./den;
    reps=reps+1;
end
end